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Abstract 

Compositional fluctuations of the binary alloy result in the corresponding 
fluctuations of the thermal conductivity of the material. During crystal growth, 
these fluctuations can significantly modify the local temperature fields at the 
liquid-solid interface. This , in turn, will affect the morphological stability of the 
growing interface. In this work, the temperature dependence of the thermal 
conductivity of the solid phase has been included into the Mullins- Sekerka 
formalism. A significant effect on the onset of the instability of planar interface 
has been predicted. It has been found, in particular, that for binary systems with 
the segregation coefficient above unity a flat interface is always unstable. The 

V 

shape of the interface fluctuation should have a single harmonic character with a 
well defined wavelength. 
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Introduction 

The stability of the crystal growth process is a major factor determining 
the structural quality of crystals used in modem technology. In fact, under 
unfavorable conditions small perturbations of process parameters can escalate 
in time and result in polycrystalline or dendritic growth. It is thus of practical 
importance to determine the critical conditions for the onset of instability and the 
wavelength of the initiating perturbation which will result in an undesirable 
growth regime. This knowledge can then be utilized to select the optimum range 
of process parameters, such as growth rate and interface temperature gradient. 
Control of the crystal growth process can be optimized, when the process 
parameters are selected correctly. Evidently, a stability criterion should be 
known in order to make such a choice. 

The existing underlying theory can be split into a hydrodynamic stability 
theory, concerning bulk fluid dynamics, and the morphological stability theory, 
which focuses attention on the evolution of the interface shape during the 
solidification process. Since pioneering work by Mullins and Sekerka [ 1 ] in 
1963 , many refinements have been made to the morphological stability theory 
and many possible effects have been included in analyses. In particular, planar 
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freezing interface stability under convection has been put forward [ 2 ], also 
interface stability problem during melting has been addressed [ 3 ]. Geometrical 
shapes other than planar have also been treated [ 4 ], eutectic solidification was ■ 
studied [ 5 ]. Also, electric and magnetic effects on morphological stability have 
been studied [ 6 ] . The effects of anisotropy of crystal growth interface kinetic 
and surface diffusion have been analyzed [ 7 ]. An important extension to 
nonlinear theory has been attempted also [ 8 ]. The radiant heat transport effect 
[ 9 ] as well as diffusion dependence on concentration [10] have been 
considered. 

In all the above analyses, a solid phase was considered as homogeneous. 
However, as the interface composition of the liquid phase fluctuates, so does the 
bulk of the growing solid phase. This, in turn, leads to perturbations in the 

V 

temperature fields. As the fluctuation of the bulk composition tan reach deep 
into a solid phase from the interface, this can create a strong fluctuation of the 
temperature field. This field, in turn, would modify the shape of the interface. 
The positive feedback of this process can create an interface breakdown, and the 
negative one will lead to the enhancement of the stability. In this paper we 
consider a binary alloy directional solidification process as a model process. We 
include the inhomogeneity of the thermal conductivity of the solid phase into 
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the linear stability analysis. The buoyancy convection effect was not included at 
present into the analysis, so our results can only be relevant to the zero-gravity 
conditions. However, the interfacial fluid dynamics correction is not significant ■ 
in many instances, and inclusion of it into analysis does not lead to a major 
modification of the Mullins- Sekerka criterion [ 2 ]. An analytic expression for 
the onset of instability has been derived, which is essentially a modified 
Mullins-Sekerka formula. The included effect can lead to enhanced stability or 
to instability, depending on the sign of the solutal gradient at the interface as 
well as on the temperature dependence of the thermal conductivity. 
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Morphological stability theory 

We consider a unidirectional solidification process of a dilute binary alloy 
with the horizontal unperturbed planar liquid- solid interface moving upwards 
with a constant velocity v 0 . It is convenient to work in the moving reference 
frame with the plane z=0 coinciding with the interface. 

We consider thermal fields T s ^ , and solute concentration fields C LS , in 
solid and liquid, respectively. These fields can be split into unperturbed and 
perturbed parts as follows: 


T L (r,t)=T L0 (z)+t L (z)exp(ot+irox) , (!) 

V 

T s (r,t)=T so (z)+t s (z)exp(at+i(ox) (2) 

C^r ,t)=C L0 (z)+c L (z)exp(at+ico x) (3) 

C s (r,t)=C so +c s (z)exp(ot+icox) (4) 

i 


where T L0 , T so , C L0 , C. 0 are unperturbed components. 

Following a standard approach to the linear stability theory, we take only one 
harmonic component in the x ( horizontal to the interface) direction of all the 
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perturbed quantities. In addition, explicit exponential time dependence exp(at) 
of these perturbations is assumed. In general, the parameter a can be complex. 
The basic question the stability theory has to answer is , at what conditions the 
parameter a is negative, or, more precisely, its real part is negative, so that the 
perturbations will die out in time. In order to answer this question, a coupled 
system of second order differential equations with related boundary conditions 
has to be solved. These equations for heat conduction and solute diffusion are 


v(K g VT s )+v 0 ar s /azr6T g /at 

(5) 

K L V 7 r L +v 0 Ps /p L dT L /az=dT L /a 

(6) 

dv ps/p^q/ez =aya 

(7) 


where k are the thermal diffusivities of solid and liquid phases and D is the 
solute diffusion coefficient. 

Boundary conditions at the solid-liquid interface can be written in the form: 


t s =t l 

(8) 

q=kC L 

(9) 

vn=n(T M +mQ-T M rK-T L ) 

(10) 
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vnL v =(-k L VT L +k s V T s )n 
vn(C s -C L )=p L /p s DV C^n 


(ID 

( 12 ) 


where n is the unit vector normal to the interface, v is the interface velocity in 
the laboratory frame, k - segregation coefficient, L v - latent heat per unit volume 
of solid, p LS - liquid, solid density, k LS - thermal conductivities of liquid, solid, 
p, - kinetic coefficient, T M - melting point of the pure substance, m - the slope of 
the liquidus line, Y - the capillary constant, K - is the interface curvature, and C s 
is solid concentration. 

The time independent stationary solution on unperturbed temperature and solute 
fields can be written in the form: 


t lo (z) = T o + 1 “ ex P(-??5; z )_ 

T so (z) = To + [ 1 - exp(-^-z)] 

Clo(z) = Coo + ^jj^C® eX P( D pT z) 


( 13 ) 

( 14 ) 

( 15 ) 
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where G L , G s denote the temperature gradients at the interface in liquid and in 
solid, respectively. From the boundary condition (10) follows the relationship 
between these gradients: 


k/v K G i=k \ ( 16 ) 

A solute concentration in liquid far away from the interface (z — * «>) is denoted 
in Eq. 1 5 by . From the boundary condition (9) one can obtain C so = C M . In 
other wards, the solute concentration in growing solid is equal the concentration 
in the melt far away from the interface. 

At this point we can proceed to the linearized equations on perturbations. In 
addition to the temperature and solute perturbations we have to introduce also a 
perturbation of the interface shape. A single Fourier component of an arbitrary 
interface perturbation can be taken in the form: 

z sL = ^e x p(at+io)x) (17) 

We assume the following expressions for k s and k. as functions of the 
perturbation c s of the solid composition: 


KsCCs)— K s (l+yc s ), k s (c s )-k s (l+y'c s ) (18) 

These relations define quantities y, y'. 
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Substitution of Eqs. (1-4) into Eqs. (5-7) and retaining only linear perturbation 
quantities yields the following system of differential equations for the 
fluctuations: 

[K s (d 2 /dz 2 -co 2 )+v 0 d/dz-o]t^(z)=-rK s c s (z)d 2 T so (z)/dz 2 -yK s dc,/dz dT S0 /dz (19) 

[K L (d 2 /dz 2 -or)+v 0 p s /p L d/dz-o]t L (z)=0 (20) 

[D(d 2 /dz 2 -co 2 )+v 0 p s /p L d/ dz-a] c L (z)=0 (2 1 ) 

We also linearize the general boundary conditions Eqs. (8-12) and obtain 
at z-0: 

Lv5z SL /5 t = k s (a/5 z +^ L 5 2 T S o/3z 2 )-k L (5 t L / 5z+ z SL 5 2T L0 /az 2 )+y'c s k s 5T S( /5z (22) 

a ZsL /a=p : mz SL ac L( /az-z SL ar L0 /az+mc L -t L +T m r dhjds ? } (23) 

az SL /a(l-k)C=kv 0 5C LO /azz SL -(l-k)v 0 c L -p L /p s Dac L /az (24) 

We utilize now the following explicit expressions (see Eqs. (13-15)): 

a 

dT L0 /dz=G L exp(-v 0 p s / ( ^z/K L ) (25) 

dC L0 /dz=G c exp(-v 0 (Vp L z/D) (26) 

dT sr /dz=G;£xp(-v 0 z/K s ) (27) 

c s (z)=kc L (0)exp(ct/v 0 z) (28) 
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where c s (z) is obtained from the expression on c L (r,t) and Eq. (9) assuming no 
diffusion of solute in the solid phase. 

Utilizing (25)-(28) , the boundary conditions can then be rewritten in the 
following form: 

(G s -G l ) 8+ t(0K(0)=0 (29) 

L v a8=k s (dt s (0)/dz-8v ( /K s G s )-k L (dt L (0)/dz-Sv () /K L G L p s /p L )+G s y , k s c s (0) (30) 

o8=(i[m8G c -8G L +nic L (0)-t L (0)-T m r8co 2 ] (31) 

a8( 1 -k)C L0 =kv 0 8G c -( 1 -k)v 0 c L (0)-p L /p s Ddc L (0)/dz (32) 

The solution to the equations (19,20,21) which is finite at z — » «> can be 
given in the form: 

c L (z)=Aexp(-pz) (33) 

. t L (z)=Bexp(-a L z) (34) 

t s (z)=Rexp(a s z)+^Aexp(X.z) (35) 

where 

P .»[i + yr^] ««> 

<37) 
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^ XykG S K S 

K3(A,2-g)2 )+ V q X-<j 

y g v 0 

A “ v 0 K S 


(38) 

(39) 

(40) 


and where 


co = a 


v o ps 
dpl 


(41) 


We substitute Eqs.(33, 34,35) into Eqs. (29-32) and , after elimination of A,B,R,8 


obtain the following condition: 


L v o+ ky Kg VqG g“ k L /K L v„G^Pg/p L -k L c(^(G +k s a s )(mG cr a/ p-G s -T m r o> 2 )= 

[G s yk\ s +k s ^(X-a s )+(k L a L +k^a s )m][a(l-k)C L0 -kv o G c ]/[p L /p s Dp-(l-k)v 0 ] 

(42) 

We restrict our analysis to stationary modes only, for which the parameter a is 
real. In case of a>0 the perturbation grows exponentially in time, and, for o<0 
there is a stable growth condition as the fluctuation dies out. The condition a=0 
will then divide the parameters space into two regions, stable, and unstable. 
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We put o=0 into Eq.(42) and obtain the stability criterion in a form given by 
Mullin and Sekerka, but with a correction term present: 


( fl +4a 2 - 1 + 2kXG(a) + Ba 2 ) - mG c ( Jl +4a 2 - 1 ) + 2kG c corr = 0 
where: 

G(a) = [G s (v 0 k s /Ks +k s as) + GL(kLaL-vok L /KLPs/pL)]/(kLaL+ksas) 


(43) 

(44) 


and 


B = T m n 


V QPS V 

dp L y 


(45) 


In the limit k s ^»D the function G(a) =(G^ s +G L k L )/(k s +k L ) which is the 
quantity commonly used in the stability criterion. 

The correction term which takes into account fluctuations of composition in solid 
has the form: 


corr = G s kk s [/ + y(i^a) 2 (l +-^-a s )]/(k L a L +k s a s ) 


(46) 


DISCUSSION 


We will examine the functional behavior of G L (a) with all other parameters 
fixed. The area above this function represents the stable condition. 
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We introduce nondimens ional functions: 

a s = K S /v 0 a s 
a L = K S /v 0 a L 
a(a) = V 1 +4a 2 - 1 
f(a)= 1-^^ + as+aL 

eps(a)= 2k 2 G c K S /vo(/ + y(-i^| 3 -) 2 (l +a s )) 


(47) 

(48) 

(49) 

(50) 

(51) 


Utilizing Eq.(l 6) we can eliminate G s and obtain from Eq.(43): 

- - . L V V 0 

-mG e ot(a)(aL + k s /k L a s ) + -r — [eps(a) + (a(a) + 2kXl + a s )]+ 
+<a(a) + 2k)Ba 2 (a L + ks/k L a s ) + G L |eps(a) + f(a)(a(a) + 2k)l = 0 


The above equation defines a function G L (a), which will have a pole when the 
denominator eps(a)+f(a)(a(a)+2k) becomes zero. Let us take a closer look at 
this term. In all practically interesting cases f(a) is always positive as is a(a). 
The function eps(a) expresses the effect of the temperature dependence of the 
thermal conductivity on the interface stability. In case of eps(a) > 0 we have a 
stabilizing correction. More interesting is the case of eps(a) < 0. It will be 
fulfilled when y,7 f >0 and G c <0 or y,y f <0 and G c >0. In such cases, eps(a) 
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tends to -oo asa-» 0. It is evident that [eps(a)+f(a)(a(a)+2k)] will cross 
zero at some small value of a. Consequently, G L (a) will have a pole at this value 
and no extremum in the G L (a) will exist. In conclusion, for eps(a) < 0 the 
planar interface is always unstable, according to the present analysis. 

For semiconducting or insulating pseudobinary alloys the slope of the 
thermal conductivity as well as the diffusivity is always negative, as the heat 
conduction in these materials is predominantly by phonons, and the contribution 
of the alloy scattering to the phonon scattering processes is significant (see, for 
example [ 11 ] ). For this class of materials the analyzed effect has a beneficial 
action of stabilization at longer wavelengths of the interface perturbation in the 
case where the distribution coefficient k<l and is destabilized at k>l. Many 
examples of such systems can readily be found. It is more difficult to find an 
example of the system for which y,y’ >0 and k<l . It appears that some of metal 
alloys of the elements with similar valencies can serve as an example. Below, 
we estimate the y, y' for a number of metal alloys. For metals with similar 
valencies the thermal conductivity of the alloy is approximately a linear function 
of the composition. Also, the density of the composition can well be 
approximated with a linear dependence. The specific heat of the metallic 
element at high temperature follows the Dulong and Petit's law, and the specific 
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heat of corresponding alloys can be obtained using the Kopp-Neumann law. We 
then obtain for the specific heat of the metallic alloy: 


Cy(x) = X * Cvo/Al + ( 1 — x) * C vq/A2 


(53) 


where C vo is the constant, about 6 cal g-atom' 1 K' 1 ; A1 , A2 are the atomic 
weights of the alloy elements , and x is the alloy composition. We also use 
similar linear relations for the density and the thermal conductivity of the alloy: 



Ps(x) = x * pi + ( 1 - x) * p2 
k s (x) = x*kl +(1 -x)*k2 


(54) 

(55) 


where the thermal conductivities and the densities kl,pl,k2,p2 correspond to the 
constituent elements. Employing the relationship between the thermal diffusivity 
and the thermal conductivity of the material: 



(56) 


we then obtain for small x: 

/ = (kl-k2)/k 2 
y = y 7 - pl/p2 - A2/A1 +2 


(57) 

(58) 


These formulas yield the following values for specific alloys: 
Sn x Bi (1 _ x) : y=6.924, y f -7.434, k=0.0023 


In x Ga ( i_ x > : 

y=2.066, 

y — 1.91 1, 
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k=0.00001 

SbBi(i_ X ). 

y=1.678, 

y -2.081, 

k=4.47 

Bi*Sb ( i. x) : 

y=-0.715, 

y--0.675, 

k-0.26 

Sn x Pb (lx ). 

y=0.503, 

y-0.893, 

k=0.59 

Cd x Zn ( i. X ) : 

y=0.041, 

y--0.165, 

k-0.18 


We also included the corresponding approximate values for the distribution 
coefficient of the alloys. The presented approach is semi-quantitative only, 
especially considering the fact that most of the above alloys form eutectics. 

As an illustration of the presented stability analysis, at Fig.l we have 
pknted the temperature gradient G L in the melt at the interface as a function of 
the wavelength of the perturbation. Also, the well known Mullin-Sekerka result 
(dotted curve) is plotted at this graph. The following set of parameters typical 
for metallic alloys has been used for this purpose: 


D=10' 5 cm 2 /sec 


T TMO- 5 cm K 

Hi } 

1^=300 joule/cm 3 
m=-3 


k=0.1 


k L =0.15 joule/(cm K sec) 
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=0.3 joule/(cm K sec) 
k l =0.1 cm 2 /sec 
k s =0.2 cm 2 /sec 

C. -0.1 

Y-Y-5 

P s /Pl =1 

v 0 =5*10' 4 cm/sec 

The pole position of G L (a) is at the wavelength of approximately 0.7 cm. 
Finally, it has been found by varying numerical values of the parameters, that 
the position of this pole is a sensitive function of the ratio D/v 0 . The smaller the 
ratio, the smaller the corresponding wavelength of the pole. Practically, this 
instability mode can potentially be detected in case of small diffusion 
coefficients and high growth rates. Also, this kind of instability, if present, will 
lead to a well defined periodic interface shape. 
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CONCLUSION 

In this work, we have investigated theoretically the effect of temperature 

$ 

dependence of the thermal conductivity on the onset of the planar interface 
instability during the crystal growth. The Mullins- Sekerka results have been 
modified and the additional correction was included into the analysis. In most of 
the cases with the distribution coefficient k smaller than one the considered 
effect yields only a small stabilizing correction to the Mullins- Sekerka formula. 
For k>l, however, our analysis shows that there is a wavelength at which the 
perturbation of the interface will always grow in time. A similar effect can be 
expected when the electric current is utilized as a driving force for the crystal 
growth (electrolytic deposition, LPEE). 
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Figure caption 


Fig.1 . The threshold temperature gradient Gl in the melt at the 
interface as a function of the wavelength of the perturbation. Solid line 
- this analysis; dotted line - Mullin-Sekerka result. 



